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Figure 1: Deformation using meshfree C 2 -weighting, (a) Handle-driven deformation based on linear blending is an intuitive method for the 
interactive shape manipulation, (b) Artifacts caused by highly distorted triangles can be generated from the weights computed on an originally 
well-meshed domain, (c) We propose a meshfree framework to generate C 2 -continuous weights for linear blending based deformation. Our 
approach inherits the merits of mesh-dependent weighting schemes meanwhile bringing the weighting method to the resolution of infinity. 


Abstract 

Handle-driven deformation based on linear blending is widely used 
in many applications because of its merits in intuitiveness, effi¬ 
ciency and easiness of implementation. We provide a meshfree 
method to compute the smooth weights of linear blending for shape 
deformation. The C 2 -continuity of weighting is guaranteed by the 
carefully formulated basis functions, with which the computation 
of weights is in a closed-form. Criteria to ensure the quality of 
deformation are preserved by the basis functions after decompos¬ 
ing the shape domain according to the Voronoi diagram of handles. 
The cost of inserting a new handle is only the time to evaluate the 
distances from the new handle to all sample points in the space 
of deformation. Moreover, a virtual handle insertion algorithm has 
been developed to allow users freely placing handles while preserv¬ 
ing the criteria on weights. Experimental examples for real-time 
2D/3D deformations are shown to demonstrate the effectiveness of 
this method. 

CR Categories: 1.3.5 [Computer Graphics]: Computational Ge¬ 
ometry and Object Modeling—Geometric algorithms, languages, 
and systems 

Keywords: shape deformation, meshfree, closed-form formula¬ 
tion, linear blend skinning 

1 Introduction 

Shape deformation techniques have various applications in com¬ 
puter graphics for image manipulation, geometric modeling and 
animation. Compared with other deformation strategies, handle- 
driven methods outperform others as they are intuitive, effective 
and easy-to-implement in many different scenarios. Using handles, 
users can bind a shape U with the handles and then manipulate their 


locations and orientations to drive the deformation of U. Specif¬ 
ically, each handle Hi with i = 1,..., m is defined as a local 
frame with its origin h 4 G O. After defining an affine transfor¬ 
mation T i for each handle Hi , the deformation of U is realized by 
computing the new position of each point p G O via a linear blend¬ 
ing of affine transformations T^p. The linear blending is weighted 
by fields Wi : U 9ft associated with handles Hi. Basically, 
to achieve an intuitive and high-quality deformation, the following 
criteria on the weights are demanded: smoothness, non-negativity, 
partition-of-unity, locality/sparsity, and no-local-maxima (see the 
analysis given in [Jacobson et al. 2011]). 

The recent advancement of technology focuses on computing 
weights of blending on a discrete form of domain (i.e., meshes are 
employed to determine piecewise linear fields of weights). Weights 
are computed on the mesh nodes via minimizing some discrete 
differential energies (e.g., biharmonic, triharmonic and quatrahar- 
monic used in [Jacobson et al. 2012b]). After incorporating the hard 
constraints according to the criteria on weights, the weights are de¬ 
termined on mesh nodes with the help of non-linear optimization. 
However, this is time-consuming. As a result, the insertion of new 
handles cannot be realized in real-time as new routines of non-linear 
optimization need to be taken. Moreover, the determined weights 
are mesh-dependent. For a symmetric shape to be deformed that is 
asymmetrically meshed, the computed weights for a handle located 
at the symmetric positions can rarely be symmetric. For poorly 
meshed computational domains, the artificial distortion caused by 
the elements of poor shape is more serious (as illustrated in Fig.l). 
Although the artifacts can be reduced by increasing the density 
of meshes, this will further slow down the computation. Ideally, 
the distribution of weights should only be affected by the shape to 
be deformed and the locations of handles, which indicates mesh- 
independence. Existing mesh-independent approaches in literature 
for handle-driven deformation (e.g., [Singh and Fiume 1998; Mil- 
liron et al. 2002; von Funck et al. 2006; Sumner et al. 2007]) can 
only satisfy subsets of the demanded properties on weights. This 
motivates our work on investigating a new meshfree method to de¬ 
termine weights for shape deformation. 

In this paper, we formulate the evaluation of weights in a closed- 
form so that the deformation framework based on this gains the 
benefit of flexibility - i.e., the response of inserting new handles 
is real-time. Specifically, the time cost of inserting a new handle 
is linear to the number of samples used to represent the domain of 
























computation. The basis function formulated in this approach can 
guarantee the properties of smoothness, non-negativity, partition- 
of-unity, locality/sparsity, and no-local-maxima, all of which are 
necessary to ensure a deformation of high-quality. 

The main results of our work are as follows: 

• We present a meshfree method to determine linear blending 
weights with C 2 -continuity for real-time deformation. The 
weights are formulated in a closed-form of basis functions 
centered at the handles (details are given in Section 3.1). Af¬ 
ter decomposing the region to be deformed by the Voronoi 
diagram of handles, aforementioned criteria of shape defor¬ 
mation are all ensured (see the analysis in Section 3.2). 

• A virtual handle insertion algorithm is proposed in Section 4 
to guarantee the locality and sparsity of weighting so that a 
deformation interpolates the transformations defined on han¬ 
dles. The virtual handles are added to let the supporting region 
of the basis function defined on a handle not cover the origins 
of any other handles (see the algorithm in Section 4.1). 

• After constructing the Voronoi diagram of all handles (includ¬ 
ing user-input and virtual ones), its dual-graph gives a con¬ 
nectivity of the handles. We compute harmonic fields on the 
graph to determine the transformations of virtual handles ac¬ 
cording to the transformations specified on the user-input han¬ 
dles (see Section 4.2). It is found that the transformations 
determined in this way lead to a shape-aware deformation fol¬ 
lowing the intention of user input. 

With the help of a discrete implementation on point samples intro¬ 
duced in Section 5, an efficient and effective meshfree approach 
has been developed for handle-driven shape deformation. 2D/3D 
experimental results are shown in Section 6 to demonstrate the per¬ 
formance of our approach. 

2 Related Work 

Shape deformation is an important research area in image manipu¬ 
lation and geometric modeling. There are a large amount of exist¬ 
ing approaches in literature. The purpose of this section is not for 
a comprehensive review. We only focus on discussing the handle- 
driven deformation approaches. 

Mesh-based techniques for discrete geometry modeling and pro¬ 
cessing have been widely explored in the past decade. Typical 
approaches including variational surface deformation [Botsch and 
Kobbelt 2004], Poisson deformation [Yu et al. 2004], Laplacian 
editing [Sorkine et al. 2004] and other linear variational surface de¬ 
formation approaches (see also the survey in [Botsch and Sorkine 
2008]). Volumetric information and rigidity are also incorporated 
to enhance the shape-preservation in [Igarashi et al. 2005; Botsch 
et al. 2006; Botsch et al. 2007; Sorkine and Alexa 2007]. One com¬ 
mon drawback of these approaches is that the positions of vertices 
on a model need to be determined by solving a system of linear 
equations after every update of handles, which becomes a bottle¬ 
neck of computation. A recent development in [Jacobson et al. 
2011; Jacobson and Sorkine 2011; Jacobson et al. 2012b] trans¬ 
fers the workload from online optimization to offline. Specifically, 
the weights corresponding to handles are computed on every vertex 
of a model before manipulating the handles (similar to [Zayer et al. 
2005]). The deformed shape is then evaluated by linear blending 
of transformations defined on handles. In [Sumner et al. 2007], the 
handles are elements of a simplified mesh. Although this strategy is 
more efficient than the deformation methods based on online opti¬ 
mization, they still cannot avoid solving large linear systems, which 
slows down the response of deformation after inserting new han¬ 


dles. Moreover, the results of deformation are also suffered from 
the artificial distortions caused by the problems of meshes (e.g., 
too coarse meshes for a fine deformation, a mesh with ‘needle’ and 
‘cap’ triangles, and the problem of symmetry). Our meshfree ap¬ 
proach solves these problems by providing closed-form formulas 
to generate weights preserving all the demanded properties for pro¬ 
ducing deformations with high quality in real-time. 

Another thread of researches for deformation focuses on mesh- 
independent approaches. Different handles are employed for shape 
manipulation. Points are used in [Yoshizawa et al. 2002; Schae¬ 
fer et al. 2006], and curves are employed as handles in [Lazarus 
et al. 1994; Singh and Fiume 1998]. Grid-based deformation tech¬ 
niques in [Sederberg and Parry 1986; Lee et al. 1995] conduct the 
bivariate/trivariate cubic splines to realize deformations with C 2 - 
continuity. Users are allowed to move control points of the spline 
surfaces/solids to modify the embedded shapes, where the editing is 
indirect. Some approaches have been developed to extend this ap¬ 
proach to provide the ability of direct editing (ref. [Hsu et al. 1992; 
Hu et al. 2001]). However, the computational domain is still limited 
to a simple topology (i.e., genus zero). An improvement of the grid- 
based techniques is introduced by Beier and Neely [1992] to allow 
handles in the form of line segments by using the Shepard’s interpo¬ 
lation [Shepard 1968]. Cage-based deformation (e.g., [Joshi et al. 
2007; Ben-Chen et al. 2009]) can be considered as a further gen¬ 
eralization of grid-based deformation, where weights can be found 
by a closed-form in terms of the handles in [Ju et al. 2005; Lipman 
et al. 2008]. However, the construction of cages is usually not au¬ 
tomatic and the manipulation on cages instead of a model itself is 
indirect. 

Moving least square (MLS) strategy is employed in [Schaefer et al. 
2006] for interpolating the similarity/rigid deformation at handle 
points. A closed-form solution is provided in their approach to de¬ 
termine the transformation matrix on every point in a MLS man¬ 
ner. The transformations in the whole domain need to be com¬ 
puted when any handle is moved. In other words, the deformation 
is globally affected by all handles - lack of sparsity. Different from 
this MLS approach, our approach belongs to the category of linear 
blending based deformation. When the property of sparsity is pre¬ 
served on the weights, the deformation at a point is only affected 
by the nearby handles that is easier to be predicted by end-users. 
Moreover, the deformation determined by our approach is resolu¬ 
tion independent, which is very important for image manipulation. 

The work of generating weights for linear blending also relates to 
the research of scattered data interpolation, where radial basis func¬ 
tions (RBF) are widely used (e.g., [Floater and Iske 1996; Botsch 
and Kobbelt 2005]). In [Botsch and Kobbelt 2005], the deformation 
is governed by global RBFs that lead to a dense linear system to be 
solved. The weights determined by the dense (or global) data inter¬ 
polation approaches lack of sparsity. Therefore, every point in the 
domain is changed when any handle is updated even if it is far away. 
Although the compactly supported radial basis functions (CSRBF) 
can help on introducing the sparsity (ref. [Floater and Iske 1996]), 
it does not provide closed-form formulas as our approach. 

3 Meshfree Weighting 

Following the linear blending formulation, the new position of a 
point p G 9 is determined by the transformations T* defined on 
handles Hi as 1 

m 

p' = Y Wi ( p ) TiP (!) 

i=1 


1 T i is a homogenous matrix and p is represented by homogeneous co¬ 
ordinate. 



( 2 ) 



Figure 2: Voronoi diagram based method to determine the size of 
local support, (a) The Voronoi diagram of handles can decompose 
Vi into smaller pieces, (b) The illustration ofruif^i) and r d (hi) in 
the Voronoi diagram, (c) Very close handles can lead to r d (hi) < 
rd(hi). 


with Wi(-) being the scalar field of weights to be determined. The 
origin of a handle Hi is denoted by h*. This linear blending based 
deformation is fast and easy-to-implement. However, carelessly as¬ 
signed weights can lead to visible artifacts in results. Basically, a 
deformation with high quality must have the following properties: 

• Smoothness: The scalar field of weights must be smooth to 
avoid visual artifact (discontinuity) in both 2D and 3D de¬ 
formations. We use compactly supported Bezier basis func¬ 
tions in our formulation, which lead to a weight field with 
C 2 -continuity. 

• Interpolation: The final transformation determined by the 
linear blending must interpolate the transformations at the 
handles. Specifically, the weight on a handle Hi is one at 
its origin while basis functions centered at other handles give 
zero at this point. This is guaranteed by the locality and the 
sparsity in our formulation. 

• Consistency: When applying the same transformation T on 
all handles, all points in ii must be consistently transformed 
by T. This is enforced by the partition-of-unity property in 
our formulation. Another consistency requirement is about di¬ 
rection. The region influenced by a handle should not change 
in the inverse direction of the transformation assigned on the 
handle. We ensure this by the property of non-negativity. 

• Shape-awareness: This is a property more or less subjec¬ 
tive. Basically, the intrinsic requirement on shape-awareness 
is to have deformations like stretching, bending and twisting 
an elastic solid, where the handles serve as pins. In our formu¬ 
lation, this is preserved by 1) having non-positive first deriva¬ 
tive of basis functions and 2) letting all basis functions have 
similar support sizes. No-local-maxima on weights will pre¬ 
vent generating singularity (e.g., a point moves faster than all 
its neighbors) during deformation. 

Our formulation below leads to C 2 -continuous weights preserving 
all these properties in deformations. 


the weights for Hi is then defined as 


m( p) 


4>i(d{p,hi)/n) 


which enforces the partition-of-unity. 

To be shape-aware and interpolate handles, fi(-) is chosen as a 
monotonically decreasing function with fi (0) = 1 and fi (t) = 0 
(Vt > 1). A quintic polynomial is employed for the function fi(t) 
so that the constraints for C 1 and C 2 -continuity at the boundary of 
the supporting regions can be satisfied. Specifically, we need 

4>'i( 0) = ^(1) = 0"(O) = ^'(1) = 0. (3) 


To ease the evaluation and analysis, each <fi ( t ) is represented as the 
^-component (i.e., fi{f) — b y (t),t G [0,1]) of a 2D Bezier curve 
with degree-n (n > 5) 


= (4) 

i =0 

where Bi, n (t) are the Bernstein polynomials. From the property 
of Bezier curves (ref. [Farin 2002]), we know that x — t when 
bf — i/n. Letting b>Q 1 2 = 1 and n -i n -2 — 0 can satisfy 
these constraints at the endpoints (see Appendix B for more details). 
For the rest control points, we can simply assign them as 0.5 or 
align them along the line b 2 b n _2 uniformly. 

When the intrinsic-distance is used to generate the input param¬ 
eter t for the basis functions, linear blending based deformations 
driven by these basis functions behave in a shape-aware manner. 
Now the problem left is how to determine the support size n of 
each basis function. As a basic requirement of handle-driven de¬ 
formation based on linear blending, every point p G O should 
be influenced by at least one handle. To be shape-aware, a point 
p should be mostly affected by its closest handle in Vi. Voronoi 
diagram sited at the origins of handles {h^} provides an intrinsic 
decomposition of Q according to these observations (see Fig.2(a)), 
where the intrinsic-distance in Q is used as the metric for generat¬ 
ing the Voronoi diagram. We denote the cell that corresponds to h ? 
by V(hi). Two metrics according to a handle Hi can be defined as 
follows (see Fig.2(b) for an illustration): 

• The size of a Voronoi cell: r d (hi) = sup qGV ( h .) d( q, hi); 

• The separation to other sites: r^(hi) = infh^ d(hi : h j). 

To let the basis function <fi{f) centered at Hi cover all points in 
V(hi) and to ensure the handle interpolation property, it should 
have 

r d (hi) <n < ^(hi). (5) 

The support size can be n = (1 — a)r d (hi) + ar^(h^) with 
a G (0,1] being specified by users as a shape factor. For most 
of the examples in this paper, a = 1 is used. It is possible to have 
two handles too close to each other so that r d (hi) < r d (hi) (see 
Fig.2(c) for an example). For solving such cases, we will use the 
virtual handle insertion algorithm (presented in Section 4). 

3.2 Analysis and Discussion 


3.1 Formulation 

Each handle Hi is equipped with a compactly supported basis func¬ 
tion with support size n as <f>i(d( p, h i)/ri), where h* is the loca¬ 
tion of Hi and d(*, •) returns the intrinsic-distance (see Appendix A 
for the definition) between two points inside Q. The scalar field of 


We analyze the advantages of our formulation for the handle-driven 
deformation based on linear blending. 

Non-negativity: 4>i(t) > 0 so that Vp G Q,Wi(p) > 0. More¬ 
over, when ri > r d (hi) is ensured for all handles, every point in Q 
should be covered by at least one handle’s support. In other words, 
Er=i^(-)#o. 


















Partition-of-unity: This has been enforced by the formulation in 
Eq.(2). That is, 


y>(p) = 

i= 1 i= 1 2 -— = 1 


(j)i(d(p,hi)/n 


'Ll (t>j(d(p,hj)/rj 


= 1 . 


Locality/Sparsity: This is preserved by Vt > 1 : fii(t) = 0 and 
the condition given in Eq.(5). The transformation at a point co¬ 
incident with a handle is only determined by the handle itself. 

Vi ^ jAji^i) = 0 . 

Smoothness: C 2 -continuity is preserved on the weights deter¬ 
mined by Eq.(2). First of all, the basis function fii(t) — b y (f) is 
C n -continuous for t E (0,1) when b y (t) is defined as a Bezier 
curve in Eq.(4) with n > 5. Therefore, Wi(p) is also C n - 
continuous when </)j(d(p,hj)/rj) ^ 0 for any other j i. In 
the region that is only covered by the support of Hi, Wi = 1. Sim¬ 
ilarly, it is also a constant function (wi = 0 ) in the region outside 
the support of Hi. By Eq.(3), it is not difficult to prove the C 2 - 
continuity at the following two cases: 

i) d( p, h j) < rj and d( p, hi) = n, 

ii) d( p, h j) m r d and d( p, hi) < n, 

where both the first and second derivatives are zero. 



Figure 3 : By using different basis functions in our formulation, 
different distributions of weights can be obtained which lead to the 
change of deformation behaviors. Isocurves for the weight field of 
the right handle are also shown in black lines on the deformed bars. 


- i.e., the influence of handles does not decay while increasing the 
distance to the handle’s center. Both the problems will be solved 
by applying the virtual handle insertion approach presented in the 
following section. 

4 Virtual Handle Insertion 


No-local-maxima: The global maxima of a weight Wi only hap¬ 
pens at the origin of handle Hi and the regions only covered by 
the support of Hi. Besides, we also observe the phenomenon of 
no-local-maxima in all our experimental tests. 

Closed-form: The weights {777 (p)} at any point p E O are evalu¬ 
ated in a closed-form (i.e., by Eq.(2)). This guarantees the flexibil¬ 
ity of inserting new handles during the deformation in real-time. 

Meshfree: As the evaluation of basis functions to determine the 
weights is only related to the intrinsic-distance from points to the 
origin of handles, the solution is independent of mesh quality and 
resolution. In the mesh-dependent solutions, elements with poor 
shape, which can occur after a drastic deformation step, must be 
optimized. Remeshing leads to another round of weights computa¬ 
tion that could be time-consuming. 

In short, our method preserves all the merits of prior methods for 
linear blending based deformation (e.g., [Jacobson et al. 2011; Ja¬ 
cobson and Sorkine 2011; Jacobson et al. 2012b]) while introducing 
new benefits of flexibility and efficiency. 

Besides the flexibility of inserting new handles during the defor¬ 
mation, we also provide users a method to change the behavior 
of handles by adjusting the shape of basis functions (i.e., fiiifi)). 
For example, as shown in Fig.3, for the basis function built 
by a septic Bezier curve (n = 7), we can assign different values 
to bf and b\ to obtain different shapes for fiiif) to have differ¬ 
ent deformation behaviors. Basically, a ‘flat’ basis function (e.g., 
= bf = 0.5) results in a deformation simulating hard materials 
while a more curved basis function (e.g., bf = 1 , bj = 0 ) makes 
the deformation soft. When using polynomials in higher orders, we 
have more degree-of-freedoms to change the shape of basis func¬ 
tion. However, according to our experiments, septic polynomials 
are good enough in most of the cases. 

The formulation of meshfree weighting also has some limitations. 
First, the interpolation property cannot be preserved when the dis¬ 
tance between two handles are too close while the regions to be 
covered by either handle are large. Specifically, the interpolation 
of handles becomes an approximation when 772 (hi) < r^(h^) in 
Eq.(5) can NOT be satisfied. Second, for the region that is only cov¬ 
ered by one handle, the transformation is consistent with the handle. 
Then, the deformation presented in this region is not shape-aware 


A handle insertion algorithm is developed to enrich our meshfree 
weighting framework in the aspects of guaranteeing the handle in¬ 
terpolation property and improving the shape-awareness of defor¬ 
mation. 

4.1 Insertion algorithm 

When 772 (hi) > 77i(hi), we know that there are points in the 
voronoi cell V(hi) whose distances to hi are larger than the mini¬ 
mal distance from hi to other handles. 

Proposition 1 When 772 (hi) > 77* (hj), inserting new sites at the 
points h d E V(hi) with d(hd,hi) = 772 (hi) can reduce fd(hi) 
while keeping 7*7 (hi) unchanged. 

Proof. First of all, the value of r^(hi) is not affected. When hd 
is the only point in V(hi) with d(hd, hi) = 772 (hi), it is obvious 
3q E V(hi) with d(q, h d ) < d(q, hi). Define S(h d ) = {q E 
V(hi) | d( q, hd) < d(q, hi)}. After inserting a new site at hd, 
the points in S(h d ) become the member of V(hd). When all points 
with d{ hd, hi) = 772 (hi) have been assigned to other voronoi cells, 
the value of 772 (hi) reduces. On the other aspect, the distances from 
the newly inserted points to hi are 772 (hi) which is greater than 

r h ( hi). □ 

Based on this proposition, we develop a greedy algorithm for han¬ 
dle insertion. Define TL as the set of handles and £(•) = 772 (•) — 
r/ l (-). When 3hi E TL with 5(hi) > 0, new handles are inserted to 
resolve this problem by reducing maxh i G'H{^(hi)/rd(hi)}. The 
pseudo-code is described as Algorithm Virtual Handle Insertion. 

Remarks. From Proposition 1, we know that inserting new han¬ 
dles in a voronoi cell V(hi) with 5(hi) > 0 can reduce the value 
of 5 (hi). However, inserting a new site hd E V(hi) can also affect 
the other handles (i.e., Hj with j 1 ). In extreme cases, the origi¬ 
nal S(hj) < 0 could be turned into S(hj) > 0. Then, new handles 
need to be added into V(h j). 

Our virtual handle insertion algorithm can be considered as a vari¬ 
ant of the farthest point sampling algorithm, which tends to tessel- 
late a domain into a voronoi diagram with neighboring voronoi cells 
having similar sizes. The condition of r d (•) < 771 (•) is satisfied on 











Algorithm 1: Virtual Handle Insertion 

Input: the set % of real handles 

Output: the expanded set Ft with virtual handles 

while eFL, 5(hQ > 0 do 

Find the handle h m = argmaxh^K 5(hi)/rd(hi); 
Find a point p £ V(h m ) with d( p, h m ) = rd(h m ); 
Insert a new handle located at p into FL; 

Update the values of r*d(-) and r^(-) on all handles; 

end 

return FL ; 



Figure 4: When two handles are too close to each other (see left), 
the condition for interpolation (i.e., rd(-) < can only be 

satisfied after inserting virtual handles. (Middle) The newly in¬ 
serted virtual handles (in blue) tessellate the deformation domain 
into voronoi cells whose areas are similar to neighboring cells. The 
Delaunay graph, VQ(FL U Ti v ), of the voronoi diagram is also 
shown - see the network linking the handles. (Right) The score, 
max{(5(hi)/rd(h^)}, of our Virtual Handle Insertion algorithm 
drops while inserting virtual handles. 


all handles when this is the case. Our experimental tests also follow 
this observation (see Fig.4 for an example). 



Figure 5: A handle covering a large region can affect the inter¬ 
polation on its nearby handles, (a) For the handle at the right, its 
voronoi cell covers all the right part of the alligator - this leads to a 
value ofrd(-) that is much larger than r^(-). In this case, transfor¬ 
mations at the left two handles cannot be interpolated, (b) Virtual 
handles (in blue color) are added to resolve the problem by the in¬ 
sertion algorithm. As a result, the domain to be deformed has been 
decomposed into smaller voronoi cells with handles (real and vir¬ 
tual) as sites, (c) The deformation result is driven by both the real 
and the virtual handles, where the transformations at real handles 
are interpolated. 


4.2 Transformation on Virtual Handles 


A left problem is how to determine the transformation on virtual 
handles according to the user-specified transformations on real han¬ 
dles. Denote the set of real handles as Fi and the set of virtual han¬ 
dles as FL V . As aforementioned, the handles of Ft U FL V have parti¬ 
tioned the given domain D into a voronoi diagram Vor(FL U Fi v ). 
A dual graph of Vor(FL U FT) can be constructed by 1) using the 
sites of every voronoi cell as nodes and 2) linking the sites of every 
two neighboring voronoi cells by a straight line, which is a Delau¬ 
nay graph [Berg et al. 2008]. We denote the Delaunay graph by 
VQ(FL U Fi v ) and also use symbol H to represent nodes in VQ 
since each node is in fact a handle (real or virtual). The transforma¬ 
tions of handles in FL V are determined with the help of the Delaunay 
graph as follows. 

• For each handle Hi in Ft , a harmonic field vo ;(•) is com¬ 
puted on VQ to assign each handle H g a field value vji(H g ). 
Boundary conditions, vji(Hi) — 1 and vji(Hj^i) = 0, are 
given to compute the harmonic field vj *(•). If there are m 
handles in FL,m harmonic fields are determined on VQ. 


• After converting the transformation T i of each handle into a 
rotation quaternion q* and a translation vector t i, the rota¬ 
tion and the translation on a virtual handle H v £ Fi v can be 
determined by 


1 


sum (H v ) 


^ ^ '&?i(H v 


Hien 

with w aU m(-) = T. Hi en w o{-)- 


( 6 ) 


• Finally, the quaternion and the translation determined on each 
virtual handle are converted back into a transformation matrix 
to be used in linear blending. 

The transformation of virtual handles determined in this way brings 
in the effect of shape-awareness during the deformation. As illus¬ 
trated in Figs.5 and 6, the deformation of whole domain driven by 
the transformations on handles (real and virtual) is very natural. 
The influence of a real handle decays when the distance to it in¬ 
creases. 

5 Implementation Details 

Similar to many other meshfree approaches, we sample the input 
domain D to be deformed into a set of dense points V. By search¬ 
ing fc-nearest-neighbors of each point, a graph Q(V) spanning D 
(in discrete form) can be established by using points in V as nodes 
and adding links between neighboring points. Note that user spec¬ 
ified handles should also be added into V to construct the graph 
(i.e., FL C V). The intrinsic-distance from any point q £ V to a 
handle is approximated by the distance between q and the handle 
on the graph, which can be computed efficiently with the help of 
Dijkstra’s algorithm. Also, the voronoi diagram Vol(FL ) can be ob¬ 
tained by the Dijkstra’s algorithm with multiple sources on Q(V ), 
where each sample is assigned to a voronoi cell. As the primitives 
used in the computation are points, the deformation approach can 
be easily generalized from 2D images to 3D solids. More examples 
can be found in the following section. To determine the weights on 
a general point p £ Q that is not a sample in V, a linear blending 
based on reciprocal distance weights [Floater and Reimers 2001] is 
























Figure 6: The deformation of a rabbit is drive by four real handles 
(see the yellow dots and the frames shown in the top row). The 
result of deformation is determined with the help of virtual handles 
(shown in blue dots). The transformations at handles (both real and 
virtual ones) are illustrated by frames. 


employed to obtain the weight on p from its /c-nearest-neighbors 
in V. There are more sophisticated parameterization strategies in 
[Floater and Reimers 2001], which can also be applied here. With 
the help of this meshless parameterization, we can easily take an 
up-sampling step in the domain TL when the point set V becomes 
sparse when applying a drastic deformation. 

After using the virtual handle insertion algorithm to generate a set 
of new handles, harmonic fields are computed on a dual graph of 
Vol(TL) to determine the transformations on virtual handles. By our 
boundary condition, all field values are non-negative when uniform 
Laplacian is employed [Wardetzky et al. 2007]. In other words, 
the coefficients used in Eq.(6) are non-negative. Instead of solving 
a linear system to compute the harmonic field, we initially assign 
the field values on all real handles as one and the weights on all 
virtual handles are set as zero. Then we apply Laplacian operators 
to update their field values iteratively. The field values on virtual 
handles can be efficiently obtained after tens of iterations. 

The point handles can be generalized to different types of handles 
(e.g., line segments and polygons, etc.). Specifically, each handle 
H g now becomes a set of points {h^} instead of a single point 
while all these points are equipped with the same transformation 
T g . The major change is the method to evaluate the intrinsic- 
distance from a query point q to handles (e.g., line segments), 
which is the intrinsic-distance to q’s closest sample point on the 
handle. The rest of our approach will keep unchanged. Extreme 
case occurs when two line-segment handles have a common end¬ 
point so that rh(-) of these two handles becomes zero. There is no 
way to satisfy the condition of raf) < r/ l (*) for handle interpola¬ 
tion. We therefore only approximate the transformations specified 
on handles. Specifically, the basis function is changed to a global 
Gaussian 

Mt) = e~ (citf (7) 

with a being a constant to control the width of Gaussian. In our 



Figure 7: A portrait is edited by segment handles. Salient feature 
inside the closed loop of segment handles at the left eye is preserved 
after the deformation. 


Figure 8: An example of processing the photograph of Tibet palace 
by segment handles. The image can be easily warped with the help 
of segment handles provided in our framework. 


implementation, letting a be ^rh(-) works well in all tests. As 
some handles may have common endpoints, r^(-) is changed to the 
minimal non-zero distance to other handles to exclude those con¬ 
nected handles. It is clear that the transformation at the position of 
a handle Hi is commonly determined by all handles in Ti although 
the influence of far away handles is trivial. On the other aspect, 
the smoothness of deformation is improved to C°°. Cages can be 
formed by linking the segment handles into closed loops. For ex¬ 
ample when editing the portrait shown in Fig.7, the cage located 
at the boundary help resize the image. Moreover, the cage at the 
left eye fully controls the shape inside it and therefore preserves the 
salient feature. 

6 Results 

Our meshfree weighting method provides a compact tool to assign 
continuous weights for all points in the domain of deformation. 
With the help of sophisticated techniques for assigning transfor¬ 
mations on the handles (e.g., the pseudo-edge method in [Jacobson 
et al. 2011] or the optimization method in [Jacobson et al. 2012a]), 
a natural user interface for shape deformation can be achieved. 

We have tested this approach in a variety of examples by using both 
the point and the segment handles. Figures 1, 5 and 6 have al¬ 
ready demonstrated the functionality of point handles. Especially, 
in Fig.5, the scheme of virtual handles insertion guarantees the in¬ 
terpolation at real handles. Figure 6 illustrates the effectiveness of 
our method in determining transformations on virtual handles. The 














Figure 9: An example of symmetric deformation: when applying 
symmetric transformations on two symmetric handles to deform a 
symmetric domain, our meshfree approach guarantees to obtain a 
symmetric result. 



Figure 10: The flexible 3D Octopus can be easily manipulated by 
using the point handles. 


example of using segment handles to deform a portrait has been 
shown in Fig.7. Another example is given in Fig.8 to warp the 
the shape of palace. To obtain natural bending results, we can add 
rotations on handles by heuristic methods (e.g., the pseudo-edge 
[Jacobson et al. 2011]). Another example is to demonstrate the per¬ 
formance of our approach in a symmetric deformation. When de¬ 
forming a symmetric domain by adding symmetric transformations 
on symmetric handles, it is expected to get a symmetric result. This 
property is preserved by our formulation (see Fig.9). 

We also apply this method to deform 3D models. In these exam¬ 
ples, the 3D models are represented by polygonal mesh surfaces. 
The weights computed by our approach are used in a linear blend¬ 
ing way to determine the new positions of vertices. Note that, the 
space enclosed by a mesh surface need to be sampled into points 
with the help of voxelization technique (e.g., [Schwarz and Sei¬ 
del 2010]) in order to evaluate the discrete intrinsic-distance in the 
domain to be deformed. Point handles are used to manipulate the 
flexible Octopus in Fig. 10, where the interface of manipulation be¬ 
comes user-friendly after employing the scheme of pseudo-edges 
to determine the transformation of point handles. Linear blend¬ 
ing scheme is widely employed in the animation of skeletal mod¬ 
els (e.g., [Jacobson and Sorkine 2011; Magnenat-Thalmann et al. 
1988]). The example shown in Fig.l 1 gives the performance of our 
approach in this scenario. 3D models with very complex topology 
(e.g., the Buddha model with internal truss structrues in Fig. 12) that 
are hard to be meshed can be easily handled in our approach. When 
deformations with large rotation are applied (e.g., in Fig. 13), a pro¬ 
gressive deformation strategy can help generate satisfactory results. 

For prior mesh-based approaches, the numerical system must be 
solved once more when new handles are inserted. In our mesh¬ 
free weighting formulation, the time cost of adding new handles is 
very trivial as the weights are determined in a closed-form. Table 
1 lists the statistics of our approach on different examples. All the 
tests are conducted on a computer with Intel Core i7-3740QM CPU 
at 2.70GHz with 8GB memory, where our current implementation 
only uses a single-core. All results of deformation can be obtained 



m 

\S\ 

tvoi (sec.) 

t w (sec.) 

Gingerman 

2(8) 

155,457 

0.584 

0.054 

Alligator 

3(3) 

53,225 

0.128 

0.015 

Rabbit 

4(15) 

22,972 

0.128 

0.015 

Portrait 

11 

4,225 

0.029 

0.005 

Palace 

22 

4,225 

0.041 

0.003 

Chinese 

2 

4,076 

0.008 

0.001 

Octopus 

10 

7,485 

0.039 

0.002 

Armadillo 

17 

26,002 

0.100 

0.014 

Buddha 

4 

236,661 

0.302 

0.051 


20 

236,661 

0.834 

0.160 

Bar 

2 

4,765 

0.009 

0.002 


Table 1: Computational Statistics for the examples shown in the 
paper. \TL\ denotes the number of handles (the number of virtual 
handles is shown in the bracket) and |«S| represents the number of 
sample points used in the computation. The columns under tvoi 
and t w state the time used in the computation of the voronoi dia¬ 
gram and the weights respectively. 


at an interactive speed. 

Discussion. When using the meshfree formulation presented in the 
paper to deform real 2D/3D objects, sample points are adopted as 
the medium for realizing the computation. The error-bound of com¬ 
putation on this discrete representation is guaranteed by the density 
of samples. However, during the process of a sequence of defor¬ 
mations, the density of points could be changed dramatically. In 
this sense, a dynamic up-sampling step should be integrated in the 
framework to preserve the error-bound of intrinsic-distance compu¬ 
tation. The image editing applications can be implemented by us¬ 
ing either the super-sampling technique or the texture mapping on 
a mesh. In our framework, the cost of weight evaluation is trivial 
after resampling. The bottleneck is the computation of intrinsic- 
distances on the sample points. Our current implementation is 
based on the Dijkstra’s algorithm. However, this shortest path prob¬ 
lem with multiple sources can be computed in parallel on the system 
with many-cores [Rong et al. 2011], which can result in a signifi¬ 
cant speedup and will be implemented in our future work. 

Our formulation gives global maximum at the positions of handles, 
which is very important to avoid the unintuitive behavior of de¬ 
formations. For a shape-aware deformation, it is also demanded 
having no-local-maximum. This has been verified in our experi¬ 
mental tests. We check the topology of isocurves on the fields of 
weights (see Fig. 14 for an example). If there is a closed loop formed 
by isocurves of Wi(-) at one place except the center of the handle 
hi(-), a local maximum is generated there. However, no such case 
is found in all our examples. 

7 Conclusion 

We present a method to determine weights of blending for shape 
deformation. Our formulation is meshfree and in a closed-form, 
which can be easily used in a variety of applications in 2D/3D de¬ 
formations. Equipped with a virtual handle insertion algorithm, 
good properties of weights generated by prior mesh-based meth¬ 
ods can all be preserved in this approach. A variety of examples 
have been shown to demonstrate the function of our approach. 

Only linear blending deformations are tested in the paper. We plan 
to further extend the application of weights generated in this ap¬ 
proach to more advanced skinning methods, such as dual quater¬ 
nion [Kavan et al. 2008], with which the blending of two rigid mo¬ 
tions will result in a rigid motion. This is a very important prop¬ 
erty when the deformation of articulated characters is computed by 
















Figure 11: An example of using the weights determined by our approach in the animation of Armadillo driven by the point and the segment 
handles. 



Figure 12: An example of deforming a model with very complex topology - the Buddha model with interior truss structures, where our 
meshfree approach can determine the weights for linear blending effectively and efficiently. 


the skinning methods. The deformations driven by linear blend¬ 
ing are not always injective and therefore can generate the results 
with foldovers and self-intersection. Recently, some researches 
have been conducted in this direction to produce injective map¬ 
pings (e.g., [Aigerman and Lipman 2013; Schuller et al. 2013]), 
which are mainly mesh-based. In a function based formulation, the 
injectivity of a mapping can be checked by the sign of Jacobian. 
However, it is still not clear about how to resolve the problem when 
self-intersection is detected. This will be one of our future work. 
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Appendix A: Intrinsic-Distance 

For a 2-manifold shape in 2D/3D Euclidean space, all points on the 
shape form a bounded domain Cl. For any two points {p s , p e } E 
Cl, if there exists a curve line C C Cl connecting p s and p e , we 
define the intrinsic-distance of {p s , p e } along the curve C as 

d(p s ,p e ]C) = length(C ) 

Then the intrinsic-distance of {p s , Pe} in Cl is defined as 
d( p s ,p e ) = rmn length{C). 

If there is no curve connecting p s and p e , that is the case they are 
not located in a connected region of Cl. The intrinsic-distance is 
then defined as d( p s , p e ) = oo. 

Sampling based intrinsic-distance. For a set of sampling points 
S G Cl of Cl, we can build a graph Q by using the sample points as 
nodes. We represent the shortest distance between p s and p e on Q 


as dg( p s , p e ;<S). If for any two points {p s , p e } G Cl, we always 
have 

|ds(Ps,p e ;<S) - d(p s ,Pe)| < e , 

the sampling S is a distance-bounded sampling of Cl. 

The intrinsic-distance defined in this way has the following proper¬ 
ties: 

• Existence: d(p s , p e ; C) is always calculable once C is deter¬ 
mined, which is a curve segment in Cl. Therefore, d(p s , p e ) 
always exists for Cl when p s and p e are located in the same 
connected region. 

• Uniqueness: d( p s , p e ) is uniquely determined while the cor¬ 
responding curves may be multiple. 

• Convergency: For any s > 0, there always exists an infinite 
sampling of Cl - that is the sampling density D(S) oo. 
Since lim £) ( 5 )^ 00 5 = 0, we have 

lim |d 0 (p s ,p e ;«S) - d(p s ,p e )| = 0. 

D(«S)—>• oo 


Appendix B: Endpoint Constraints 


From the analysis in [Farin 2002], we know that 

b'(0) n(bi - b 0 ), b'(l) n( b n - b n -i) 

for a Bezier curve in n-th order. And also 


b /7 (0) = n(n — l)(b 2 — 2bi + bo) 
b"(l) » n(n - l)(b n - 2b n _i + b n _ 2 ) 
Incorporating the constraints in Eq.(3), we have 


bi b 0 ,b n b n _i,bi 


bo + b 2 
2 


, b n _i 


b n + b n _ 2 
2 


As we already need bf = i/n to let x = t, it is not difficult to find 
that bg = b^ = b| = 1 and b^ - b y n _i — b y n _ 2 — 0 satisfy all 
these constraints. 



